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ABSTRACT 

Using interferometric data from BIMA observations, combined with detailed 
modeling in Fourier space of the physical structures predicted by models, we 
constrain the circumstellar envelope parameters for four Class young stellar 
objects, as well as their embedded circumstellar disks. The envelopes of these 
objects are still undergoing collapse, and theoretical collapse models can be com- 
pared to the observations. Since it has been suggested in a previous study that 
both the Larson-Penston and Shu similarity solutions underestimate the age of 
the system, we adopt Tassis & Mouschovias' model of the collapse process, which 
includes all relevant magnetic fields effects. The results of the model fitting show 
a good consistency between theory and data; furthermore, no age problem exists 
since the Tassis & Mouschovias' model is age independent for the first 255 kyr. 
Although the majority of the continuum dust emission arises from the circumstel- 
lar envelopes, these objects have well known outflows, which suggest the presence 
of circumstellar disks. At the highest resolution, most of the large-scale envelope 
emission is resolved out by interferometry, but the small-scale residual emission 
remains, making it difficult to observe only the compact disk component. By 
modeling the emission of the envelope and subtracting it from the total emission, 
we constrain the disk masses in our four systems to be comparable to or smaller 
than the typical disk masses for T Tauri systems. 

Subject headings: stars: format ion — stars: pre-main-sequence — circumstellar 
matter — radio continuum: stars — techniques: interferometric — magnetic 
fields 
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Introduction 



The standard scenario of low-mass star formation starts at the collapse of prestellar 
cores and the formation of central protostellar objects. These young stellar objects (YSOs) 
evolve through the so-called Class 0, I , II, and HI stages, which are thought to be a temporal 



Lada fc Wilkind 11984 ; lAdams et al 



1987; Andre et al. 



1993 . 12000). In the 



sequence (e.g. 

earliest stage, i.e. Class stage, when the cen tral YSO is just form ing inside the surrounding 
envelope (of mass ~ a few solar masses, e.g., lLooney et al.ll2000l . hereafter LMW2000), the 
envelope is still undergoing gravitational collapse onto the circumstellar disk. The YSO 
powers the bipolar outflows, which carve away the polar regi on of the envelope by entraining 
envelope material and wide ning their opening angles (e.g., lBachillerlll996l ; lArce fc Sargent 



2006 



Seale &: Looneyj 2007 ). At this early time, the envelope mass is >85% of the system 



mass (lLooney et al.ll2003l . hereafter LMW2003). As the system evolves, the envelope loses 
mass as material is transported down through the circumstellar disk onto the protostars or 
carried away with the outflows. Eventually, the YSO circumstellar structure is dominated 
by the disk (a hundredth of a solar mass, e.g., Andrews & Williams 2005). The circumstellar 
disk evolves, presumably becoming a planetary system like the Solar System. 

The initial collapse process of low-mass protostars is o ften described by self-similar 
isothermal solutions, which are a continuum of soluti ons (e.g.. Iwhitworth fc Summers 1985 ) 
that range from the "inside-out " collapse so l ution (Shu 19771, hereafter the Shu solution) 



.ap. 



1969; Penston 



196 



9; 



Hunter 1977. hereafter the LP 



to the Larson-Penston solution (lLarsor 
solution). These models generally obtain an inner core with a power law density profile 
p oc r -3 / 2 that increases in radius with time, surrounded by a p oc r~ 2 envelope. The 
theoretical density profiles from these sol utions have been compared to observatio ns of the 
dust continuum emission (e.g., LMW2003; lHarvey et al.ll2003l ; |j0rgensen et al.ll2005l ). but the 
models cannot fit the observations with reasonable physical parameters (required age is too 
low; see LMW2003), which consequently hints at the need for more sophisticated theoretical 
models that include more of the essential physical processes of the collapse mechanisms, for 
example, turbulence and/or magnetic fields. 

The theory of turbulence- induced star formation postulates that turbulence causes over- 
densities and is thus responsible for the core formation in molecular clouds, while magnetic 
fields are not d ynamically important and do not have a significan t impact on this process 
(see reviews of iMac Low fc Klessenll2004l ; lElmegreen fc Scald 12004 ). At this time, there are 
no predictions of the density of a protostellar object that is produced by turbulence-induced 
collapse. Moreover, the non-thermal contribution of the observed linewidths is small in 



evolved, collapsing mo lecular cloud cores (e.g. jBenson fc Myerdll989l ; lBarranco fc Goodman 
19981 ; iKirk et al.l 120071 ) . Whether turbulence plays an important role in the formation and 
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evolution of protostellar fragments (or cores) is still under debate. On the other hand, 
the theory of ambipolar-diffusion-initiated star formation predicts the formation of ther- 
mally and magneticall y supercritical protostellar cores inside magnetica lly subcritical parent 
clouds (see reviews of iMouschoviad Il996l ; iMouschovias fc Ciolekl Il999l ). In the framework 
of the ambipolar-diffusion -induced collapse, there are e xtensive studies of the dynamics of 
the prestellar phase (e.g., Tassis fc Mouschovias 2007a b cj) and the accretion p rocess after 
a protostar has formed at the center of the core (ITassis fc Mouschoviasl l2005al ibl. hereafter 
TM2005), which can be compared to observations. 

TM2005 have constructed a six-fluid MHD simulation to study the accretion process of 
matter from a molecular cloud core onto a protostellar object in the presence of magnetic 
fields. In their model, they track the evolution of magnetic flux and six kinds of particles 
(neutral molecules, atomic and molecular ions, electrons, neutral grains, negatively-charged 
grains, and positively-charged grains, among which only the electrons are assumed to be 
attached to the magnetic field lines) in a self-gravitating, accreting molecular cloud core. 
The simulation starts with a magnetically supported parent cloud. Ambipolar diffusion 
leads to the formation of a thermally and magnetically supercritical core that begins to 
contract dynamically. Its innermost part reaches nearly hydrostatic equilibrium while its 
outer part still undergoes infall. At the moment when the hydrostatic protostellar core has 
just formed at the center (called t = 0), the inner core region including the protostar is 
replaced by a central sink to facilitate the calculation. As time progresses and mass and 
magnetic flux accrete onto the protostar from the envelope, a region of enhanced magnetic 
field, called the "magnetic wall", forms and drives an outward-propagating shock. Behind 
the shock, gravity dominates over the magnetic forces and reaccelerates the neutral particles, 
which continue to accrete onto the protostar until the next magnetic wall is formed. The 
magnetic wall forms and disperses in a quasi-periodic manner. Because of the presence of 
the magnetically controlled bursts, accretion from the envelope onto the protostar is episodic 
over the first 255 kyr. 

Interferometric data of Class objects provide the best means to test these theories. 
One of the features of an interferometer is the ability to spatially filter emission. Indeed, dust 
continuum observations of young protostars have been often used to peer inside the bright 
envelope to reveal the young, compact circumstellar disk (e.g-. lKeene fc Massonlll990l ). Dust 
continuum emission is often used, but molecular lines can also be excellent tracers of specific 
conditions. However, using molecular lines to trace the disks do have some difficulties in the 
case of the youngest stars: (1) active accretion and outflow processes at multiple scales, (2) 
chemistry effects and evolution, and (3) shocks in the outflow, the disk, and the disk/envelope 
interface region. All of these contaminate the preferred disk-only tracers and make it difficult 
to disentangle any molecular result without a good understanding of the envelope structure 
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derived from the dust continuum modeling (e.g.. iBrinch et al.ll2007l ). 



Regardless of the tracers used, more intricate theoretical models than the "inside-out" 
collapse can be tested observationally. In this paper, we build on the observational data of 
LMW2000: A = 2.7 mm dust continuum images of 24 young stellar sources with sensitivity 
to spatial scales from (X'5 to 50". A discussion of the data acquisition and images can be 
found in that paper. We follow the work of LMW2003, which presented modeling of the 
envelope emission of Class objects, and use the 4 sources from that work with the highest 
signal-to-noise ratio (NGC 1333 IRAS 4A, NGC 1333 IRAS 4B, NGC 1333 IRAS 2A, and 
L1448 IRS 3B). We compare the predicted density profiles from TM2005 to these data and 
comment on the results with respect to those found in LMW2003. 



Sources 



All 4 sources modeled in this study are in the Perseus molecular cloud, a low mass 
star forming region probably in the vicinity of the massive star forming region Per OB2 
association. Dense cores and YSO candidates at all evolutionary stages (Class 0, I, II, and 
III) have been found in Perseus via radio and infrared observations. IC 348 and NGC 1333 
are the two main dense clusters and other smaller groups like L1448, L1455, Barnard 1 , 
and Barnard 5 are al so associated with many low mass protostars (e.g. Enoch et al.l 12006 



Jprgensen et al.l 120061 ) . The exact distance to the Perseus molecular cloud is still uncertain 
and ranges from 220 to 350 pc. Th e smaller val ue is based on the distance-interstellar 
extinction relation using photometry (jCernisI Il990l ); the distance may be the larger value 
if the Perseus molecular cloud is p hysically related to th e Per OB2 association, Hipparcos 
parallax distance of 318 ± 27 pc (jde Zeeuw et al.lll999l ). Since Perseus is composed of a 
long chain of dense clouds with a total length of about 30 pc, t here may be a distance 
gradient or it may be composed of several layers of clouds. ICernisI (119931 ) showed a distance 
difference from the eastern part (~260 pc) to the western part (~220 pc). In this study, we 
adopt a distance of 350 pc for NGC 1333 and 300 pc for L1448 as in LMW2003 to facilitate 
comparison. Since we are using roughly the upper limit of the distance, this assumption may 
lead to overestimates of the envelope and disk masses and underestimates of the source size. 



2.1. NGC 1333 IRAS 4 



NGC 1333 is a reflection nebula with mainly two embedded protostellar clusters in the 
L1450 dark nebula in the Perseus molecular cloud. The age of the young cluster is about 1 
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Myr, estimat ed by the fraction of infrared excess sources and a K-band luminosity function 
com parison (ILada et al.lll996f) . consistent with the age estimated by brown dwarf studies 



.e.g. 



Wilking et al.ll2004l ). The plentiful jet and outflow activities driven by YSOs in this 
region also imply that it is an active sta r forming region at an early stage of evolution (e.g., 



Bally et al.lll996l : Knee & Sandelll l2000h . 



The multiple system NGC 1333 IRA S 4 (hereafter IRAS 4) contains mainly three groups 
of sources designated as 4A, 4B, and 4C (ISandell fc Kneel 1200 ll ). The brightest Class com- 
ponent IRAS 4A has b een detected as a binary system separated by 1.8" with a common 
circumbinary envelope (jLay et al.lll995t LMW2000) A highly collimated N-S molecular out- 
flo w driven by IRAS 4A2 has b een observed with estimated dynamical age of about 6000 
yr (IBlake et al.l Il995l ; I Choi! 120051 ) . A dimmer southern outflow probably driven by 4A1 has 



been mapped in HCN and SiO, but no northern counterpart has been detected (jChoil 12001 



20051 ). SMA polarimetric observations have shown the magnetic field geo metry of IRAS 4A 
and supported the magnetic theory of star formation (IGirart et al.l 120061 ) . 



IRAS 4B is a binary system with a separation of ~10" between 4BW and 4BE (ISandell fc Knee 



200ll ); 4BE is also named as 4B', 4BII, or 4C in different references. Note that 4C is mostly 
used as the name of another millimeter object ~50" n ortheast of 4A. A compact collimated 
outflow driven by IRAS 4BW has been observed (e.g., 
and shown a short dynamical timescale. 



Choill200ll ; |Pi Francesco et al.l 1200 if ) 



2.2. NGC 1333 IRAS 2A 



NGC 1333 IRAS 2 (hereafter IRAS 2) has been resolved into 3 sources ( ISandell fc Knee 



2001) including two Class protostars, IRAS 2A and IRAS 2B, and a starless core IRAS 
2C (I J0rgensen et al.l 120041 ) . Two CO outflows hav e been mapped: one in the NNE-SSW 



direction and the other in the E-W direction (e.g., iKnee fc Sandelll 120001 ) . Since these two 
outflows are orthogonal to each other and have quite different properties, they may have 
different driving sources. The E-W o utflow may be driven b y IRAS 2A, while the NNE-SSW 
outflow may be driven by IRAS 2C (IKnee fc Sandelll 120001 ) . It's also possible that both are 
driven by IRAS 2A which may be an unresolved binary. IRAS 2B may drive the third ou tflow 



in this region but only the blue-shifted lobe has been identified (IKnee fc Sandelll |2000| ). 
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2.3. L1448 IRS 3B 



The star forming region L1448 is located ~ 1° southwest of NGC 1333 in Perseus and 
contains many YSOs: Class I source L1448 IRS 1, Class sources L1448 IRS 2, L1448 IRS 
3 (al so known as L1448 N) , and L1448-mm (also know as L1448 C) are the most well-known 



(e.g. Uorgensen et al.ll2006l ; lO'Linger et al.ll2006l ). Many of these YSOs have been resolved 
into binary systems; for example, L1448 IRS 3 and L1448-mm are both binary systems. 
Multiple outflows have been found in this region. The huge molecular outflows emanating 



Wolf- Chase et al. 


2000; 


Kwon et al. 


2006 



Tobin et all 120071 ). 



L1448 IRS 3 is composed of 3 sources, among which IRS 3A (L1448 N:A) and IRS 3B 
(L1448 N:B) are separated by 7" and have a common envelope in a protobinary system, 
and IRS 3C (L1448 NW) is 20" northwest of them (LMW2000). All of them are Class 
objects, except that IRS 3A is slightly closer to the transition between Class and Class I 
(b'Linger et alJbood ). Two interacting o utflows driven by IRS 3A and IRS 3B, respectively, 
have been studied by lKwon et al. ( 2006 ). 



3. Modeling 

To compare the theoretical models of TM2005 and the LMW2000 observations, we 
characterize the observed dust continuum emission, which depends on the dust density, 
grain properties, and temperature. 



3.1. Density Profiles 

In the theoretical TM2005 model, the physical parameters of the envelope are shown 
to have a quasi-periodic variation beginning at the time of about 15 kyr after the formation 
of a hydrostatic protostellar core. The density goes through a cycle profile that is largely 
invariant with time from 15 to 255 kyr. This implies that, unlike the LP or Shu models, we 
cannot estimate the age of the source based on its density profile alone. 

In this study, we adopt a typical set of density profiles, averaged along the characteristic 
scale height, from one of the magnetic cycles predicted by the TM2005 model. The density 
structure repeatedly goes through phase = 1 to 15 in an evolution cycle, as shown by the 
solid curves in Fig. [TJ Here we only plot the more representative phases. The dotted curve 
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in Fig. [T] shows the initial density profile at the time hydrostatic equilibrium is established 
in the central region of the core (t = 0), which nearly follows a power-law relation of index 
— 1.7. We can choose density profiles of any cycle for this study, as they are all similar. 
In other words, different cycles have the same predicted density profiles and quantitatively 
similar phases. The age of the chosen magnetic cycle is from t = 33,750 to 37,250 yr after 
the formation of the central hydrostatic core. By this time, the series of magnetic cycles has 
been well established. 

The model cloud in TM2005 has an equatorial radius of 4.23 pc. However, the thermally 
and magnetically supercritical fragment extends to approximately 9000 AU and contains ~9 
M . Outside this region, hardly any evolution takes place (see Fig. [1]). In order to compare 
with observations, we truncate the density profiles at an outer radius within which the 
dynamical infall takes place. This is the "envelope", at the center of which the protostar 
is forming. By scaling the density profiles, the theoretical TM2005 model is applied to 
observed objects of different masses. (This kind of scaling of the density profiles implies a 
corresponding scaling of the magnetic field profiles, but, since no information on the magnetic 
field is available from observations, we do not discuss this implication further in this paper.) 



Although the TM2005 model has a flattened morphology (IFiedler &: Mouschoviaslll993l ). 
we only use a spherical density profile. The flattened molecular cloud has a radius of several 
pc, which is much larger than the scale studied here. LMW2003 has tested geometric effects 
by elongating the envelope, and shown that flattening the envelope will artificially make a 
steeper flux density profile in u-v space. However, the analysis normalized the maximum 
flux density, which is different than the effect studied here. To better compare the effect 
a flattened envelope has on our modeling, we modified our spherical envelope model by 
multiplying the predicted image with a flat-shaped mask with exponentially-decayed edges 
similar to a flattened envelope. The half-thickness of the "disk" is one-fifth of the radius, 
which gives a similar flattening ratio to that in the TM2005 model. In this case, the observed 
flux density is less steep at small u-v distances, but is unchanged beyond u-v distances of 
10 kX. This is because fewer large-scale components, or small u-v distance components, 
contribute in the flattened structure compared to a spherical model. We fit the simulated 
data of a flattened envelope with our spherical model and found that the spherical model 
fits the data very well with the same density profile, but that the mass is overestimated. In 
addition, the density profile from TM2005 has been averaged along the characteristic scale 
height, which also helps mitigate the effect of a flattened envelope. On the other hand, the 
vertical density pro file needs to be better modeled (see the observed flattened envelope in 



Looney et al.l 120071 ). which is beyond the scope of this paper. Finally, a circumstellar disk 
may exist inside the envelope and its emission is also taken into account in the model; this 
is described in more detail in 5 3.4 and § 4 . 
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3.2. Grain Properties and Temperature Profiles 



Properties of dust grains such as composition and size determine their extinction and 
emission, as well as the temperature stru cture and the observable d ust emission of the en- 
velope. Here we adopt the grain model of IWolfire fc Cassinelfil (119861). in which the material 



mixtu re and size distributions of the Mathis-Rumpl-Nordsieck model (IMathis. Rumpl. fe Nordsieck 
1977 . hereafter called the MRN model) and optical constants of grains in Draine fc Led ( 1984 ) 
are used. The grain model consists a mixture of uncoated graphite and silicate with particle 
sizes ranging from 0.005 \im to 0.25 fim and a power law size distribution of index -3.5. It 
is important to note that the grains in TM2005 are spherical and uniform in size with a dif- 
ferent chemical composition from MRN. However, these grain d ifferences are not significant 
in the theoretical model evolution (IDesch fc Mouschovias!l200ll ). 



The mass opacity coefficient k v of the grain is frequency dependent and typically follows 
a power law relation k v oc . The index varies with environment and is related to 
grain properties. At submillimeter wavelengths, the unevolved grains of the ISM have ~ 
2. However, in disks and dense cores de creases to 1 mainly due to grain growth (e.g. 
Beckwith fc Sargentlll99ll ; iNatta et al.ll2007l ). Although our observations are only at a single 
wavelength, k u of a wider range of frequency is still needed for computing a self-consistent 
temperature profile since radiation of all frequencies contributes to the total luminosity. Our 
model uses the MRN grain model = 2 at optical and infrared wavelengths but assumes 
— 1 at long wavelengths, and adopts k v = 0.009 cm 2 g^ 1 at A = 2.7 mm. The model 
also assumes the dust grain properties are uniform and do not change with radius in the 
envelope. 

Temperature profiles are then considered based on the model grain properties. The tem- 
perature profile can be simplified as a power law with radius if the dust envelope is optically 
thin and the dust opacity has a power-law frequency dependence {k v oc u 13 ): temperature 
T oc r~ 2 '^ + ^' assuming the central protostar is the only heating source. In the case of — 1, 
the temperature is T oc r~ ' 4 . But in our model, the inner part becomes optically thick as 
the density increases near the center (LMW20 03), so we calculate a self-c onsistent temper- 
ature profile for each fitting using the code of IWolfire fe Cassinellil (119861 ). The code takes 
the luminosity of the central object and solves the radiative transfer equation including the 
effects of both emission and extinction by the dust grains in the envelope. At each shell of 
the envelope, the luminosity is conserved. 

The self-consistent temperature profiles are important for calculating the emission, but 
they contradict the isothermal assumptions in theoretical models. However, until non- 
isothermal theoretical models are developed, this is the best compromise (also see LMW2003). 
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3.3. Interferometric Filtering 

Although we use the power of the interferometer to resolve out the large-scale features of 
the envelope, it is important to point out that there is remnant envelope emission even with 
high resolution configurations (e.g., LMW2000). There are a few reasons for this. First, 
the envelope emission is p ower-law like (e.g., LM W2003), so the expected interferometric 



response from a Gaussian (IWilner fc Welch! Il994f ) is not applicable. Second, the inner edge 



cutoff is abrupt, especially with the steep density profile. The abruptness causes ripples in 
u-v space that create power at long baselines. The important point here is that the ability 
to detect a disk in the presence of an envelope is not set by the formal noise level or Fourier 
components, but by the intrinsic ability to model the complexity of the envelope (inner 
cutoff, asymmetric structure, etc.) and the resolution of the observations compared to the 
disk size and the inner cutoff of the envelope. In other words, it is difficult to separate the 
youngest circumstellar disk from the inner circumstellar envelope; revealing the embedded 
circumstellar disk still requires an understanding of the inner envelope. 

Only interferometric data are used in this study. Single dish data can give an upper 
limit at zero-spacings, but should not be included in the fitting. Single dish observations 
detect not only the emission from the inner envelope, but also that from the outer part of 
a cloud that is not actively involved in the protostar forming process. Using interferometry, 
only structures of size scale similar or smaller than the collapsing envelope are detected and 
modeled. Hence the mass and radius we are inferring in this paper are not the total mass or 
the overall size of the prestellar cores, but that of the inner envelope undergoing the collapse, 
observed at A = 2.7 mm. 



3.4. Modeling Details 

The observational data we use are from LMW2000. To compare the observations and 
the theoretical model, we analyze the data in u-v space where the data are not affected by 
the CLEAN algorithm or u-v sampling. The interferometric data are binned in u-v annuli 
around the source locations from LMW2000 and averaged vectorially. The resulting u-v 
amplitudes for each bin are shown by asterisks in Fig. [2J The displayed error bars are 
statistical error bars based on the standard deviation of the mean of the data points in the 
bin with a minimum of 10%, reflecting the uncertainty in the overall calibration. In the 
cases where the binary systems were separated by more than 10", the companion sources 
were subtracted out of the u-v data using the large-scale images of LMW2000. The new u-v 
data were remapped to confirm that the large-scale emission from the companion sources 
was not detected. Although there may exist some residual of small-scale emission from 
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the companion envelope in the u-v data, vector averaging in u-v annuli will minimize its 
contribution. 

The observational data shown in Fig. [2] show the brightness distribution plotted at 
various antennas, or u-v, spacings, which is the Fourier transform of the sky brightness 
distribution. In other words, power at small u-v distance represents large-scale structures 
and power at large u-v distance represents small-scale structures, i.e., a point source would 
be a constant at all antenna separations. The brightness distribution is determined by the 
circumstellar material surrounding the source depending on density, temperature and grain 
properties. 

We model the circumstellar envelope emission as arising from a spherically symmetric 
dust envelope with TM2005 radial density and self-consistent temperature profiles and an 
embedded circumstellar disk represented by an envelope-attenuated point source. The ob- 
served emission of most sources shows a circular symmetry and lack of significant internal 
structures (see images of the sources in LMW2000). The combined radiative transfer code 
allows the calculation of the expected flux as a function of radius in the image plane or as a 
function of u-v distance in the Fourier plane. This provides the best way to trace emission 
structure to very small length scales, effectively the density and temperature profiles in the 
inner circumstellar envelope. 

For each object, we did a parameter fitting of the model to the observational data with 
four degrees of freedom: evolutionary phase (density profiles) in the TM2005 magnetic cycle, 
outer cutoff radius of the envelope, central point source flux, and envelope mass. The inner 
cutoff radius is fixed to be 10 AU to be consistent with the central sink approximation 
used in TM2005 and also the inner envelope is truncated by the central disk physically. In 
Fig. [3] we show the evolutionary phases in u-v space with all other parameters fixed for an 
example case of outer radius 5000 AU, envelope mass 5 M Q , no point source flux, and a 
fixed power- law temperature profile of index -0.4 and 500 K at 1 AU as an optically thin 
case (the solid curves). The difference between the phases in u-v space is mainly at large 
u-v spacings, corresponding to the shock propagation in the density profiles. On the other 
hand, adding a point source is like adding a constant in u-v space, which effects the slope 
in logarithmic plots; increasing the envelope mass increases amplitude in all u-v spacings, 
especially in short u-v spacings; changing the outer radius alters the overall level and also 
the distribution of u-v amplitudes. The effects from varying each single parameter intertwine 
together and degeneracy makes it difficult to point out which parameter is the key to each 
good fit; moreover, the self-consistent temperature profile considers the increase of optical 
depth near the center of the envelope and influences the predicted emission. 

Model images of the envelope are computed with consideration of the envelope emission 
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and the central point source attenuated by the dust envelope. The self-consistent tem- 
perature profiles are calculated from the assumed luminosity, which was derived from the 
far-infrared flux density (see LMW2003) and listed in column 2 in Table [TJ The model im- 
ages are multiplied by the observational primary beam, Fourier transformed into visibilities, 
and sampled with the same u-v coverage as the observations. Both model and observational 
data are binned and averaged in u-v annuli, shown as the flux density as a function of u-v 
spacing, and compared to each other. A reduced Xr is computed to determine the goodness 
of a fit. Among the four degrees of freedom, the total envelope mass is adjusted to minimize 
Xr while the other three parameters are fixed. This nonlinear minimization is done for any 
combination of these three parameters with outer radius ranging from 2000 to 9000 AU and 
point source flux as the ranges given in the parentheses in Table HJ and for each set of pa- 
rameters a total envelope mass was obtained with local minimal Xr- Collecting the results 
of sets of parameters, the best fit, with global minimal Xr is found. 



4. Results and Discussion 

In Fig. [21 data are binned, averaged, and shown by asterisks with associated error 
bars, and the curves show the best fit for each source as a function of u-v distance with 
the best fit parameters given on the plot. The best fit parameters do not need to be the 
same as those in LMW2003 since a different theoretical model is used here. Fig. H] shows 
the corresponding density and temperature profiles for these best fits. The straight dashed 
lines in the temperature plots are lines with slope —0.4, which would be the temperature 
profiles if the envelopes are optically thin (T oc r~ a4 ). We can see that the outer part of the 
envelopes is nearly optically thin, and becomes hotter and optically thicker at smaller radii. 

We refer to a fit with more than 90% confidence level as an acceptable fit, and a summary 
of acceptable fits is given in Table [H Fitting with density profiles at different phases implies 
that the systems may be at different phases of a magnetic cycle. Presumably, if we had 
data of more objects we might catch systems at all phases. The acceptable fits spread 
over a range of parameters depending on the signal-to-noise ratio of the observations. One 
important aspect of the modeling is the clear need for high signal-to-noise data. The x 2 
value is the evaluation quantity of the goodness of a set of parameters, and is smaller when 
the observational uncertainties are larger; low signal-to-noise ratio data make acceptable fits 
easier. When the signal-to-noise ratio is too low, the fitting becomes meaningless, which is 
why we apply the analysis only to four sources from LMW2000. This argument is consistent 
with our results that the best constraints of model parameters is provided by the highest 
signal-to-noise observation IRAS 4A among all four sources. Given better signal-to-noise 
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ratio data at more wavelengths in the future we may be able to constrain the model better 
not only for more sources but also for better-constrained parameters. 

Nonetheless, the most important result of this study is that with the TM2005 theoretical 
model, we can easily fit the observations of these sources without conflicts of ages, unlike 
the fits using the LP or Shu models (LMW2003), although the exact source age is not 
determinable by comparison with the TM2005 model. Due to the episodic nature of the 
theoretical model, we are not fitting nor implying a specific age for the source, but the range 
of age is 15 to 255 kyr. Parsec-scale jets and outflows have been found in these regions and 
can be used to estimate the ages of central sources that drive the outflows. Fo r example, 



the timescale of the outflows in NGC 1333 is of order 10 5 yr ( Bally et al.lll996l ). which is 
consistent with the age range of the TM2005 model. 

We use the phases of a typical magnetic cycle for the fitting. TM2005b Fig. 5 shows 
the beginning evolution after the central sink is introduced and the system has the first 
magnetic cycle at about t = 15 kyr; then it repeats the magnetic cycles as in TM2005b Fig. 
8 and 9 until t = 255 kyr when the sink mass reaches 1 M Q and the simulation stops. Once 
the cycling behavior has been well established, all physical quantities show high similarity 
of variation from cycle to cycle. Only the period of a cycle, which is about several thousand 
years, decreases with time, controlled by the ambipolar-diffusion timescale at the position of 
the magnetic wall; other than that, there is no obvious difference between cycles. So given 
a fit to a density profile, the system may be at a specific phase of any cycle during this 
epoch. However, the age range suggested by the theoretical model is consistent with the age 
estimated by the outflows' scale, which implies that the observational data are consistent 
with the TM2005 picture of early star formation. 

One of the most important differences between the TM2005 and the Shu or LP models 
is that the steep power-law-like density profiles of TM2005 are actually not in equilibrium, 
while p oc r~ 2 in the other models corresponding to the singular isothermal sphere is a 
critical equilibrium situation. In the Shu or LP models, the collapse is induced by an 
outward-moving rarefaction wave and the density changes dramatically after being affected 
by the wave; in TM2005, the system is collapsing and mass keeps accreting onto the central 
protostar dynamically without making an abrupt change, except for the relatively small 
bump generated by the magnetic wall. Periodic creation and dispersion of the magnetic 
wall dominate the variation of density profiles and also the accretion, so the whole collapse 
process is regulated by magnetic forces. Accounting for magnetic fields in the theoretical 
models, the density profiles at later times of the evolution are very different from those of 
the Shu or LP models, and are better matched with the observations and estimated ages 
from outflows. 
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Can we constrain the earliest disks? Tabled] shows that in all cases, we can fit our sources 
without a circumstellar disk component at the 90% confidence level. On the other hand, a 
disk is expected early on in the collapse process due to the rotation and/or magnetic fields 
in the initial cloud. As material accretes, the disk receives more mass from the envelope. 
The mass is processed in the disk, which regulates mass flow through the disk and onto 
the protostar. The magnetic fie l ds wi thin the disk and star give rise to an outflow (e.g., 
Shang et al.l 120071 ; iPudritz et al.l 120071 ) that is typically seen even in very young sources, 
suggesting that a disk is established nearly concurrently with the protostar's growth in 
luminosity as it collapses to stellar size. In fact, simple theoretical argume nts suggest that 
the d isk evolves early and grows quickly with time as ~ t 2 or ~ t 3 (e.g., 
20051 ). which depends crucially on the mass accretion rate. 
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Although not a statistical requirement, we do expect these sources to have some sort 
of deeply embedded circumstell ar disks in the center of the infalling envelopes. We can use 
HL Tau, which was modeled in iMundy et al.l (119961 ). as a standard candle to estimate the 
mass of the disks in our modeled systems. Also using BIMA A = 2.7 mm dust continuum 
observations, the total flux from HL Tau was ~100 mJy, and the derived disk mass from 
modeling was ~0.05 M . Based on the assumption of equal flux to disk mass ratio corrected 
by the distance, the value 0.010 Jy in the best fit point source flux for NGC 1333 IRAS 4A 
corresponds to a disk of mass 0.03 M , and 0.035 Jy for IRAS 4B corresponds to a 0.11 M 
disk. A distance to NGC 1333 of 350 pc has been assumed here. If we use a distance of 
250 pc instead, the estimated mass of embedded disks for the best fits would be 0.016 M & 
and 0.056 M & for IRAS 4A and IRAS 4B, respectively. The acceptable fits with maximum 
disk components give disk mass of 0.16 M , 0.18 M , 0.065 M , and 0.041 M for NGC 
1333 IRAS 4A, IRAS 4B, IRAS 2A, and L1448 IRS 3B, respectively. It gives an upper limit 
of disk mass to IRAS 2A since model parameters with higher point source flux cannot fit 
anymore. 



Unlike I Jprgensen et al.l (120051 ). a disk component is not crucially necessary in our mod- 
eling. The major reason is that different envelope models are used. TM2005 predict a very 
different visibility amplitude profile from what is expected by a power-law density profile. 
For example, in Fig. [31 the dashed curve is the predicted u-v amplitude of a nearly power- 
law density profile of index -1.7 (the initial density profile in the TM2005 simulation), and 
the dotted curve is generated with the same density profile, but with a point source flux 
representing an unresolved circumstellar disk. As can be seen, the envelope emission from a 
power-law density profile shows a very different shape than those predicted by the phases of 
the magnetic cycl e. TM2005 density pr ofiles are able to fit the data well without adding a 
Gaussian disk (cf. I Jprgensen et al.ll2005l ). Again this shows the importance of understanding 
the collapsing envelope in order to understand the embedded disk. 
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5. Summary 

1. Although interferometry is a powerful tool in resolving out the large-scale emission 
of the envelope, the ability to detect a disk in the presence of an envelope is not 
set by the formal noise level or Fourier components, but by the intrinsic ability to 
model the complexity of the envelope (inner cutoff, asymmetric structure, etc.) and 
by the resolution of the observations compared to the disk size and the inner cutoff 
of the envelope. In other words, reliable detection of the embedded circumstellar disk 
requires a knowledge of the physical parameters of the inner envelope. 

2. Our observational data are consistent with the theoretical predictions of TM2005 con- 
cerning the density profiles. Moreover, there is no discrepancy in age based on the 
size of the outflows, unlike the fitting results of the simpler Shu or LP solutions to 
the isothermal sphere (e.g., LMW2003). It is important to note that the exact age 
cannot be determined by comparison of the data with the TM2005 model due to its 
periodic nature. Regardless of the collapse initiation, the magnetic fields may play 
an important role early on in the collapse process. To expand this comparison, we 
will continue to observe Class objects at higher sensitivity and multiple wavelengths, 
better incorporating the theoretical models into our comparisons. 

3. Although our initial results do not require the existence of circumstellar disks (accept- 
able fits of 0.0 to 0.11 Mo), we can place upper limits on the disk masses. In general, 
the disks are less massive than ~0.1 M Q . The youngest circumstellar disk mass is not 
overly massive compared to other well known Class II or III circumstellar disks. The 
fact that a disk component is not statistically nec essary in this modeling is different 



from similar work at sub- millimeter wavelengths by iJprgensen et al.l (120051 ) . The main 
reason is that different visibility profiles are predicted by TM2005 than for envelopes 
of power-law density. 
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radius (AU) 

Fig. 1. — The number density profiles of neutral particles in a typical magnetic cycle at differ- 
ent phases from the TM2005 model. These curves show phase 1, 2, 3, 5, 8, 11, and 14, which 
correspond to t = 33750, 34000, 34250, 34750, 35500, 36250, and 37000 yr after formation of 
the central protostar. The shock driven by the "magnetic wall" forms, propagates outward, 
and disperses throughout a cycle. Although this is a particular cycle, no obvious differences 
are found between cycles, that is, every magnetic cycle goes through similar phases, except 
that the period of the cycle varies as the system evolves (see TM2005a,b for more details). 
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Fig. 2. — The flux density of the observational data and the best fit for each source with 
best fit parameters listed. The x 2 values are 1.49, 0.22, 1.30, and 0.41 for NGC 1333 IRAS 
4A, IRAS 4B, IRAS 2A, and L1448 IRS 3B, respectively. 
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Fig. 3. — The solid curves show the example u-v visibilities from an optically thin envelope 
with outer radius 5000 AU, mass 5 M & , and various evolutionary phases (no point source flux 
is added, and observational parameters such as the distance and u-v samplings are assumed 
to be the same as NGC 1333 IRAS 4A). Different colors correspond to different phases in 
the TM2005 model, as the density profiles of these phases are plotted in Fig. [TJ The dashed 
curve is the flux density of the same case but with the initial density profile of TM2005, 
which is a nearly power-law profile with index -1.7. The dotted curve is the same nearly 
power-law envelope emission plus a constant 175 mJy (a point source flux) representing an 
unresolved disk. The same visibilities are also plotted in a logarithmic scale in the inset 
panel. The observational data points of NGC 1333 IRAS 4A are shown for comparison. The 
power-law envelope itself (the dashed curve) is not a good fit; also, just adding a point source 
flux to the power-l aw density profil e cann ot fit the data (the dotted curve). This is similar 
to the results from iJorgensen et al.l (120051 ). However, the predicted envelope emission from 
the TM2005 model is very different and a better fit to the data. 
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Fig. 4. — The density and temperature profiles of the best fit for each source. The straight 
dashed lines in the temperature plots show the optically thin case. 
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Table 1. Acceptable Fits 



Source 


Luminosity 


Density Profiles 


Outer Radius 


Point Source 


Envelope Mass 


Good Fits d 


Name 


(L©) a 


(Phase) b 


(AU) 


Flux (mjy) c 


(Ms) 




NGC 1333 IRAS 4A 


16 


9-11 


4000-5000 


0-50 (0-50) 


4.88-6.23 


15/720 


NGC 1333 IRAS 4B 


5.2 


1-14 


2000-6000 


0-56 (0-56) 


1.45-4.02 


302/1215 


NGC 1333 IRAS 2A 


30 


2-13 


5000-8000 


0-7 (0-21) 


0.76-1.37 


49/480 


L1448 IRS 3B 


6.8 


3-12 


2000-9000 


0-18 (0-18) 


0.75-4.37 


223/480 



Note. — A set of parameters gives a good fit (a so-called acceptable fit) if the reduced \ 2 is within the 90% confidence 
level. For each parameter the range of acceptable fits is given. 

a All sources here are binary system and flux of single component is assumed based on the ratio of fluxes at A = 2.7 mm 
(see LMW2003). 

b The number corresponds to a specific phase in the cycle for which the density profile of the TM2005 model provides an 
acceptable fit (see Fig. [TJ. 

c The parameter search ranges fitting are listed in the parentheses. 

d The numbers of good fits and the total number of combinations while fitting. 



